This notebook performs pair-wise comparisons of qPCR gene expression, normalized to GAPDH expression. It calculates delta Cq, delta delta Cq, and fold changes in expression. Additionally, it generates box plots (delta Cq), and bar plots (fold change expression).
Groups are named in the following fashion:
<life.stage>.<conditioning.treatment>.<acute.treatment>
This allows for parsing downstream.
NOTE: Below is the full set of groups for the entire experiment. For the current qPCR analysis, seed and spat do not have acute treatments; just conditioning treatments.
seed.control.ambient=c("29", "40", "55", "63", "69", "101", "119", "122", "155", "164", "187", "202", "209", "214", "233", "236", "275")
seed.control.high=c("42", "59", "60", "62", "86", "102", "140", "176", "177", "184", "192", "223", "234", "243", "244", "254", "264")
seed.treated.ambient=c("14", "48", "66", "72", "89", "115", "129", "138", "156", "182", "191", "201", "227", "239", "270", "277", "280")
seed.treated.high=c("15", "19", "24", "88", "92", "105", "111", "113", "120", "128", "161", "200", "211", "256", "257", "266", "285")
spat.control.ambient=c("11", "30", "36", "52", "77", "114", "134", "142", "144", "183", "193", "229", "230", "231", "240", "272", "287")
spat.control.high=c("27", "74", "93", "96", "97", "137", "143", "153", "168", "178", "189", "206", "262", "274", "282", "284", "289")
spat.treated.ambient=c("9", "13", "38", "46", "47", "121", "145", "151", "174", "194", "197", "198", "216", "235", "241", "252", "291")
spat.treated.high=c("6", "25", "50", "78", "124", "126", "131", "160", "163", "172", "220", "226", "242", "253", "296", "298")
juvenile.control.ambient=c("18", "57", "65", "75", "79", "104", "110", "123", "125", "171", "175", "205", "238", "273", "279", "293", "317")
juvenile.control.high=c("12", "39", "43", "49", "71", "130", "141", "146", "150", "170", "195", "297", "301", "324", "351", "355", "371")
juvenile.treated.ambient=c("1", "34", "64", "83", "98", "147", "152", "158", "162", "169", "188", "271", "295", "310", "357", "361", "381")
juvenile.treated.high=c("28", "53", "61", "73", "81", "106", "109", "139", "149", "173", "181", "213", "290", "302", "311", "364", "392")
adult.control.ambient=c("3", "5", "13*", "16", "17", "80", "87", "94", "148", "159", "179", "180", "250", "258", "268", "312", "326", "330", "334", "346", "360", "377", "379", "386")
adult.control.high=c("20", "23", "26", "32", "33", "67", "70", "90", "107", "132", "135", "157", "166", "186", "207", "215", "248", "316", "341", "344", "349", "382", "394", "395")
adult.treated.ambient=c("7", "31", "35", "37", "41", "54", "84", "100", "112", "116", "118", "133", "154", "199", "203", "204", "208", "219", "294", "318", "339", "353", "363", "378")
adult.treated.high=c("21", "22", "45", "82", "85", "91", "95", "99", "103", "108", "117", "127", "165", "185", "190", "196", "232", "237", "245", "263", "276", "306", "343", "374")
# Combine vectors into lists
# Used for adding treatment info and/or subsetting downstream
groups_list <- list(juvenile.control.ambient = juvenile.control.ambient,
juvenile.control.high = juvenile.control.high,
juvenile.treated.ambient = juvenile.treated.ambient,
juvenile.treated.high = juvenile.treated.high,
adult.control.ambient = adult.control.ambient,
adult.control.high = adult.control.high,
adult.treated.ambient = adult.treated.ambient,
adult.treated.high = adult.treated.high,
seed.control.ambient = seed.control.ambient,
seed.control.high = seed.control.high,
seed.treated.ambient = seed.treated.ambient,
seed.treated.high = seed.treated.high,
spat.control.ambient = spat.control.ambient,
spat.control.high = spat.control.high,
spat.treated.ambient = spat.treated.ambient,
spat.treated.high = spat.treated.high)
Normalized to designated normalizing gene
calculate_delta_Cq <- function(df) {
df <- df %>%
group_by(Sample) %>%
mutate(delta_Cq = Cq.Mean - Cq.Mean[Target == "GAPDH"]) %>%
ungroup()
return(df)
}
# Function to create box plots for each comparison
create_boxplot_delta_Cq <- function(data, comparison, t_test_results) {
# Extract life stages from comparison
life_stages <- unlist(strsplit(comparison, "\\."))
# Debugging: Print life stages
# print(paste("Life stages for comparison:", comparison))
# print(life_stages)
# Filter data for the relevant life stages
filtered_data <- data %>%
filter(life.stage %in% life_stages)
# Debugging: Print filtered data
# print("Filtered data:")
# print(filtered_data)
# Check if both life stages are included
if (!all(life_stages %in% unique(filtered_data$life.stage))) {
stop("Not all life stages are included in the filtered data")
}
y_limits <- range(filtered_data$delta_Cq, na.rm = TRUE)
# Debugging: Print y_limits
# print("Y limits:")
# print(y_limits)
# Filter t_test_results for the current comparison
t_test_results_filtered <- t_test_results %>%
filter(comparison == !!comparison)
# Debugging: Print filtered t_test_results
# print("Filtered t_test_results:")
# print(t_test_results_filtered)
# Filter t_test_results for asterisks
t_test_results_with_asterisks <- t_test_results_filtered %>%
filter(asterisk != "")
# Debugging: Print t_test_results_with_asterisks
# print("t_test_results_with_asterisks:")
# print(t_test_results_with_asterisks)
formatted_title <- paste0(toupper(substring(life_stages[1], 1, 1)), substring(life_stages[1], 2),
" vs. ",
toupper(substring(life_stages[2], 1, 1)), substring(life_stages[2], 2))
boxplot <- ggplot(filtered_data, aes(x = Target, y = delta_Cq, fill = life.stage)) +
geom_boxplot(position = position_dodge(width = 0.75)) +
theme_minimal() +
theme(legend.position = "right") +
scale_fill_manual(values=c("darkgray", "salmon", "lightblue", "lightgreen")) +
ylim(y_limits) +
labs(x = "Target", y = "Delta Cq", title = formatted_title) +
# Highlighted section: Adds asterisks
geom_text(data = t_test_results_with_asterisks,
aes(x = Target, y = y_limits[2] - 1, label = asterisk),
vjust = -0.5, size = 8, color = "black", inherit.aes = FALSE)
print(boxplot)
}
life_stage, treatment1, and treatment2).filtered_data.The t_test_results_filtered data frame is filtered for the specific comparison.
The t_test_results_with_asterisks data frame is created to include only the rows with asterisks.
The formatted_title variable is created by capitalizing the first letter of each component and concatenating them with ” - ” and ” vs. ” in between.
This should create box plots comparing conditioning treatments within each life stage, with titles formatted as <life.stage> - Treated vs. Control.
# Function to create box plots for each comparison of conditioning treatments within life stages
create_boxplot_conditioning <- function(data, comparison, t_test_results) {
# Extract life stage and conditioning treatments from comparison
comparison_parts <- unlist(strsplit(comparison, "\\."))
life_stage <- comparison_parts[1]
treatment1 <- comparison_parts[2]
treatment2 <- comparison_parts[3]
# Debugging: Print life stage and treatments
# print(paste("Life stage and treatments for comparison:", comparison))
# print(c(life_stage, treatment1, treatment2))
# Filter data for the relevant life stage and conditioning treatments
filtered_data <- data %>%
filter(life.stage == life_stage, conditioning.treatment %in% c(treatment1, treatment2))
# Debugging: Print filtered data
# print("Filtered data:")
# print(filtered_data)
# Check if both treatments are included
if (!all(c(treatment1, treatment2) %in% unique(filtered_data$conditioning.treatment))) {
stop("Not all treatments are included in the filtered data")
}
y_limits <- range(filtered_data$delta_Cq, na.rm = TRUE)
# Debugging: Print y_limits
# print("Y limits:")
# print(y_limits)
# Filter t_test_results for the current comparison
t_test_results_filtered <- t_test_results %>%
filter(comparison == !!comparison)
# Debugging: Print filtered t_test_results
# print("Filtered t_test_results:")
# print(t_test_results_filtered)
# Filter t_test_results for asterisks
t_test_results_with_asterisks <- t_test_results_filtered %>%
filter(asterisk != "")
# Debugging: Print t_test_results_with_asterisks
# print("t_test_results_with_asterisks:")
# print(t_test_results_with_asterisks)
# Format the title
formatted_title <- paste0(toupper(substring(life_stage, 1, 1)), substring(life_stage, 2),
" - ",
toupper(substring(treatment1, 1, 1)), substring(treatment1, 2),
" vs. ",
toupper(substring(treatment2, 1, 1)), substring(treatment2, 2))
boxplot <- ggplot(filtered_data, aes(x = Target, y = delta_Cq, fill = conditioning.treatment)) +
geom_boxplot(position = position_dodge(width = 0.75)) +
theme_minimal() +
theme(legend.position = "right") +
scale_fill_manual(values=c("darkgray", "salmon")) +
ylim(y_limits) +
labs(x = "Target", y = "Delta Cq", title = formatted_title) +
# Highlighted section: Adds asterisks
geom_text(data = t_test_results_with_asterisks,
aes(x = Target, y = y_limits[2] - 1, label = asterisk),
vjust = -0.5, size = 8, color = "black", inherit.aes = FALSE)
print(boxplot)
}
life_stage, treatment1, and treatment2).filtered_data data frame is filtered to include only the rows with the relevant life stage and acute treatments.filtered_data.The t_test_results_filtered data frame is filtered for the specific comparison.
The t_test_results_with_asterisks data frame is created to include only the rows with asterisks.
Format the Title:
<life.stage> - Ambient vs. High.# Function to create box plots for each comparison of acute treatments within life stages
create_boxplot_acute <- function(data, comparison, t_test_results) {
# Extract life stage and acute treatments from comparison
comparison_parts <- unlist(strsplit(comparison, "\\."))
life_stage <- comparison_parts[1]
treatment1 <- comparison_parts[2]
treatment2 <- comparison_parts[3]
# Debugging: Print life stage and treatments
# print(paste("Life stage and treatments for comparison:", comparison))
# print(c(life_stage, treatment1, treatment2))
# Filter data for the relevant life stage and acute treatments
filtered_data <- data %>%
filter(life.stage == life_stage, acute.treatment %in% c(treatment1, treatment2))
# Debugging: Print filtered data
# print("Filtered data:")
# print(filtered_data)
# Check if both treatments are included
if (!all(c(treatment1, treatment2) %in% unique(filtered_data$acute.treatment))) {
stop("Not all treatments are included in the filtered data")
}
y_limits <- range(filtered_data$delta_Cq, na.rm = TRUE)
# Debugging: Print y_limits
# print("Y limits:")
# print(y_limits)
# Filter t_test_results for the current comparison
t_test_results_filtered <- t_test_results %>%
filter(comparison == !!comparison)
# Debugging: Print filtered t_test_results
# print("Filtered t_test_results:")
# print(t_test_results_filtered)
# Filter t_test_results for asterisks
t_test_results_with_asterisks <- t_test_results_filtered %>%
filter(asterisk != "")
# Debugging: Print t_test_results_with_asterisks
# print("t_test_results_with_asterisks:")
# print(t_test_results_with_asterisks)
# Format the title
formatted_title <- paste0(toupper(substring(life_stage, 1, 1)), substring(life_stage, 2),
" - ",
toupper(substring(treatment1, 1, 1)), substring(treatment1, 2),
" vs. ",
toupper(substring(treatment2, 1, 1)), substring(treatment2, 2))
boxplot <- ggplot(filtered_data, aes(x = Target, y = delta_Cq, fill = acute.treatment)) +
geom_boxplot(position = position_dodge(width = 0.75)) +
theme_minimal() +
theme(legend.position = "right") +
scale_fill_manual(values=c("darkgray", "salmon")) +
ylim(y_limits) +
labs(x = "Target", y = "Delta Cq", title = formatted_title) +
# Highlighted section: Adds asterisks
geom_text(data = t_test_results_with_asterisks,
aes(x = Target, y = y_limits[2] - 1, label = asterisk),
vjust = -0.5, size = 8, color = "magenta", inherit.aes = FALSE)
print(boxplot)
}
# Function to create box plots for each comparison of acute treatments within life stages and conditioning treatments
create_boxplot_acute_conditioning <- function(data, comparison, t_test_results) {
# Extract life stage, conditioning treatment, and acute treatments from comparison
comparison_parts <- unlist(strsplit(comparison, "\\."))
life_stage <- comparison_parts[1]
conditioning_treatment <- comparison_parts[2]
treatment1 <- comparison_parts[3]
treatment2 <- comparison_parts[5]
# Filter data for the relevant life stage, conditioning treatment, and acute treatments
filtered_data <- data %>%
filter(life.stage == life_stage, conditioning.treatment == conditioning_treatment, acute.treatment %in% c(treatment1, treatment2))
# Check if both treatments are included
if (!all(c(treatment1, treatment2) %in% unique(filtered_data$acute.treatment))) {
stop("Not all treatments are included in the filtered data")
}
y_limits <- range(filtered_data$delta_Cq, na.rm = TRUE)
# Filter t_test_results for the current comparison
t_test_results_filtered <- t_test_results %>%
filter(comparison == !!comparison)
# Filter t_test_results for asterisks
t_test_results_with_asterisks <- t_test_results_filtered %>%
filter(asterisk != "")
# Format the title
formatted_title <- paste0(toupper(substring(life_stage, 1, 1)), substring(life_stage, 2),
" - ",
toupper(substring(conditioning_treatment, 1, 1)), substring(conditioning_treatment, 2),
" - ",
toupper(substring(treatment1, 1, 1)), substring(treatment1, 2),
" vs. ",
toupper(substring(treatment2, 1, 1)), substring(treatment2, 2))
boxplot <- ggplot(filtered_data, aes(x = Target, y = delta_Cq, fill = acute.treatment)) +
geom_boxplot(position = position_dodge(width = 0.75)) +
theme_minimal() +
theme(legend.position = "right") +
scale_fill_manual(values=c("darkgray", "salmon")) +
ylim(y_limits) +
labs(x = "Target", y = "Delta Cq", title = formatted_title) +
# Adds asterisks
geom_text(data = t_test_results_with_asterisks,
aes(x = Target, y = y_limits[2] - 1, label = asterisk),
vjust = -0.5, size = 8, color = "black", inherit.aes = FALSE)
print(boxplot)
}
# Get a list of all CSV files in the directory with the naming structure "*Cq-Results.csv"
cq_file_list <- list() # Initialize list
cq_file_list <- list.files(path = cqs_directory, pattern = "Cq-Results\\.csv$", full.names = TRUE)
# Initialize an empty list to store the data frames
data_frames_list <- list()
# Loop through each file and read it into a data frame, then add it to the list
for (file in cq_file_list) {
data <- read.csv(file, header = TRUE)
data$Sample <- as.character(data$Sample) # Convert Sample column to character type
data_frames_list[[file]] <- data
}
# Combine all data frames into a single data frame
combined_df <- bind_rows(data_frames_list, .id = "data_frame_id")
str(combined_df)
'data.frame': 2192 obs. of 17 variables:
$ data_frame_id : chr "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" ...
$ X : logi NA NA NA NA NA NA ...
$ Well : chr "A01" "A02" "A03" "A04" ...
$ Fluor : chr "SYBR" "SYBR" "SYBR" "SYBR" ...
$ Target : chr "ATPsynthase" "ATPsynthase" "ATPsynthase" "ATPsynthase" ...
$ Content : chr "Unkn-01" "Unkn-01" "Unkn-02" "Unkn-02" ...
$ Sample : chr "206" "206" "220" "220" ...
$ Biological.Set.Name : logi NA NA NA NA NA NA ...
$ Cq : num 26.7 26.7 25.8 25.9 25.1 ...
$ Cq.Mean : num 26.7 26.7 25.9 25.9 25.1 ...
$ Cq.Std..Dev : num 0.0455 0.0455 0.0239 0.0239 0.0813 ...
$ Starting.Quantity..SQ.: num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Log.Starting.Quantity : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Mean : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Std..Dev : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Set.Point : int 60 60 60 60 60 60 60 60 60 60 ...
$ Well.Note : logi NA NA NA NA NA NA ...
# Remove rows with Sample name "NTC"
combined_df <- combined_df[combined_df$Sample != "NTC", ]
# Replace values in the Target column
combined_df$Target <- gsub("Cg_GAPDH_205_F-355_R \\(SR IDs: 1172/3\\)", "GAPDH", combined_df$Target)
combined_df$Target <- gsub("Cg_ATPsynthase_F/R \\(SR IDs: 1385/6\\)", "ATPsynthase", combined_df$Target)
combined_df$Target <- gsub("Cg_cGAS \\(SR IDs: 1826/7\\)", "cGAS", combined_df$Target)
combined_df$Target <- gsub("Cg_citrate_synthase \\(SR IDs: 1383/4\\)", "citrate synthase", combined_df$Target)
combined_df$Target <- gsub("Cg_DNMT1_F \\(SR IDs: 1510/1\\)", "DNMT1", combined_df$Target)
combined_df$Target <- gsub("Cg_HSP70_F/R \\(SR IDs: 598/9\\)", "HSP70", combined_df$Target)
combined_df$Target <- gsub("Cg_Hsp90_F/R \\(SR IDs: 1532/3\\)", "HSP90", combined_df$Target)
combined_df$Target <- gsub("Cg_VIPERIN_F/R \\(SR IDs: 1828/9\\)", "viperin", combined_df$Target)
str(combined_df)
'data.frame': 2180 obs. of 17 variables:
$ data_frame_id : chr "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" ...
$ X : logi NA NA NA NA NA NA ...
$ Well : chr "A01" "A02" "A03" "A04" ...
$ Fluor : chr "SYBR" "SYBR" "SYBR" "SYBR" ...
$ Target : chr "ATPsynthase" "ATPsynthase" "ATPsynthase" "ATPsynthase" ...
$ Content : chr "Unkn-01" "Unkn-01" "Unkn-02" "Unkn-02" ...
$ Sample : chr "206" "206" "220" "220" ...
$ Biological.Set.Name : logi NA NA NA NA NA NA ...
$ Cq : num 26.7 26.7 25.8 25.9 25.1 ...
$ Cq.Mean : num 26.7 26.7 25.9 25.9 25.1 ...
$ Cq.Std..Dev : num 0.0455 0.0455 0.0239 0.0239 0.0813 ...
$ Starting.Quantity..SQ.: num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Log.Starting.Quantity : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Mean : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Std..Dev : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Set.Point : int 60 60 60 60 60 60 60 60 60 60 ...
$ Well.Note : logi NA NA NA NA NA NA ...
levels(as.factor(combined_df$Target))
[1] "ATPsynthase" "cGAS" "citrate synthase" "DNMT1"
[5] "GAPDH" "HSP70" "HSP90" "viperin"
# Filter out rows where Cq.Std..Dev is NA
combined_df <- combined_df[!is.na(combined_df$Cq.Std..Dev), ]
# Filter rows where Cq.Std..Dev is greater than 0.5
high_cq_std_dev <- combined_df[combined_df$Cq.Std..Dev > 0.5, ]
# Print the filtered rows with specified columns, without row names
print(high_cq_std_dev[, c("Target", "Sample", "Cq", "Cq.Std..Dev")], row.names = FALSE)
Target Sample Cq Cq.Std..Dev
HSP70 244 33.10339 4.0838809
HSP70 244 27.32791 4.0838809
HSP90 223 24.85319 0.7548714
HSP90 223 25.92074 0.7548714
viperin 223 30.30089 0.6058663
viperin 223 31.15772 0.6058663
viperin 243 32.57817 0.5527617
viperin 243 33.35989 0.5527617
DNMT1 296 31.21374 0.6417578
DNMT1 296 30.30616 0.6417578
DNMT1 298 35.68716 0.5406704
DNMT1 298 34.92253 0.5406704
DNMT1 223 32.24089 0.6214201
DNMT1 223 33.11971 0.6214201
DNMT1 243 36.63921 0.5125743
DNMT1 243 35.91432 0.5125743
DNMT1 285 33.63443 0.7036122
DNMT1 285 34.62949 0.7036122
GAPDH 316 23.94926 8.5684728
GAPDH 316 24.14183 8.5684728
GAPDH 316 38.88564 8.5684728
GAPDH 213 26.98012 2.2910353
GAPDH 213 23.00009 2.2910353
GAPDH 213 26.95634 2.2910353
GAPDH 263 22.42154 0.8731474
GAPDH 263 23.77008 0.8731474
GAPDH 263 24.05667 0.8731474
citrate synthase 230 24.44066 4.4783429
citrate synthase 230 24.40421 4.4783429
citrate synthase 230 32.17909 4.4783429
viperin 227 30.47773 3.5152533
viperin 227 30.37738 3.5152533
viperin 227 36.51553 3.5152533
viperin 245 26.05748 5.1635899
viperin 245 34.98192 5.1635899
viperin 245 26.01928 5.1635899
viperin 341 26.48675 2.9838590
viperin 341 31.67235 2.9838590
viperin 341 26.52174 2.9838590
viperin 344 29.98184 2.3712440
viperin 344 25.90358 2.3712440
viperin 344 25.84648 2.3712440
viperin 355 28.79712 0.5821437
viperin 355 29.57428 0.5821437
viperin 355 28.43490 0.5821437
# Group by Sample and Target, then filter out the outlier replicate
combined.fitered_df<- combined_df %>%
group_by(Sample, Target) %>%
filter(abs(Cq - mean(Cq, na.rm = TRUE)) <= Cq.Std..Dev)
# Print the filtered data frame
str(combined.fitered_df)
gropd_df [1,520 × 17] (S3: grouped_df/tbl_df/tbl/data.frame)
$ data_frame_id : chr [1:1520] "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" "lifestage_carryover/data/qPCR/Cq/sam_2024-03-25_06-10-54_Connect-Quantification-Cq-Results.csv" ...
$ X : logi [1:1520] NA NA NA NA NA NA ...
$ Well : chr [1:1520] "A01" "A02" "A03" "A04" ...
$ Fluor : chr [1:1520] "SYBR" "SYBR" "SYBR" "SYBR" ...
$ Target : chr [1:1520] "ATPsynthase" "ATPsynthase" "ATPsynthase" "ATPsynthase" ...
$ Content : chr [1:1520] "Unkn-01" "Unkn-01" "Unkn-02" "Unkn-02" ...
$ Sample : chr [1:1520] "206" "206" "220" "220" ...
$ Biological.Set.Name : logi [1:1520] NA NA NA NA NA NA ...
$ Cq : num [1:1520] 26.7 26.7 25.8 25.9 25.1 ...
$ Cq.Mean : num [1:1520] 26.7 26.7 25.9 25.9 25.1 ...
$ Cq.Std..Dev : num [1:1520] 0.0455 0.0455 0.0239 0.0239 0.0813 ...
$ Starting.Quantity..SQ.: num [1:1520] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Log.Starting.Quantity : num [1:1520] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Mean : num [1:1520] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ SQ.Std..Dev : num [1:1520] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
$ Set.Point : int [1:1520] 60 60 60 60 60 60 60 60 60 60 ...
$ Well.Note : logi [1:1520] NA NA NA NA NA NA ...
- attr(*, "groups")= tibble [760 × 3] (S3: tbl_df/tbl/data.frame)
..$ Sample: chr [1:760] "200" "200" "200" "200" ...
..$ Target: chr [1:760] "ATPsynthase" "DNMT1" "GAPDH" "HSP70" ...
..$ .rows : list<int> [1:760]
.. ..$ : int [1:2] 21 22
.. ..$ : int [1:2] 245 246
.. ..$ : int [1:2] 85 86
.. ..$ : int [1:2] 53 54
.. ..$ : int [1:2] 117 118
.. ..$ : int [1:2] 149 150
.. ..$ : int [1:2] 213 214
.. ..$ : int [1:2] 181 182
.. ..$ : int [1:2] 321 322
.. ..$ : int [1:2] 889 890
.. ..$ : int [1:2] 637 638
.. ..$ : int [1:2] 1047 1048
.. ..$ : int [1:2] 1205 1206
.. ..$ : int [1:2] 479 480
.. ..$ : int [1:2] 731 732
.. ..$ : int [1:2] 1363 1364
.. ..$ : int [1:2] 323 324
.. ..$ : int [1:2] 891 892
.. ..$ : int [1:2] 639 640
.. ..$ : int [1:2] 1049 1050
.. ..$ : int [1:2] 1207 1208
.. ..$ : int [1:2] 481 482
.. ..$ : int [1:2] 733 734
.. ..$ : int [1:2] 1365 1366
.. ..$ : int [1:2] 325 326
.. ..$ : int [1:2] 893 894
.. ..$ : int [1:2] 641 642
.. ..$ : int [1:2] 1051 1052
.. ..$ : int [1:2] 1209 1210
.. ..$ : int [1:2] 483 484
.. ..$ : int [1:2] 735 736
.. ..$ : int [1:2] 1367 1368
.. ..$ : int [1:2] 327 328
.. ..$ : int [1:2] 895 896
.. ..$ : int [1:2] 643 644
.. ..$ : int [1:2] 1053 1054
.. ..$ : int [1:2] 1211 1212
.. ..$ : int [1:2] 485 486
.. ..$ : int [1:2] 737 738
.. ..$ : int [1:2] 1369 1370
.. ..$ : int [1:2] 329 330
.. ..$ : int [1:2] 897 898
.. ..$ : int [1:2] 645 646
.. ..$ : int [1:2] 1055 1056
.. ..$ : int [1:2] 1213 1214
.. ..$ : int [1:2] 487 488
.. ..$ : int [1:2] 739 740
.. ..$ : int [1:2] 1371 1372
.. ..$ : int [1:2] 1 2
.. ..$ : int [1:2] 225 226
.. ..$ : int [1:2] 65 66
.. ..$ : int [1:2] 33 34
.. ..$ : int [1:2] 97 98
.. ..$ : int [1:2] 129 130
.. ..$ : int [1:2] 193 194
.. ..$ : int [1:2] 161 162
.. ..$ : int [1:2] 331 332
.. ..$ : int [1:2] 899 900
.. ..$ : int [1:2] 647 648
.. ..$ : int [1:2] 1057 1058
.. ..$ : int [1:2] 1215 1216
.. ..$ : int [1:2] 489 490
.. ..$ : int [1:2] 741 742
.. ..$ : int [1:2] 1373 1374
.. ..$ : int [1:2] 333 334
.. ..$ : int [1:2] 901 902
.. ..$ : int [1:2] 649 650
.. ..$ : int [1:2] 1059 1060
.. ..$ : int [1:2] 1217 1218
.. ..$ : int [1:2] 491 492
.. ..$ : int [1:2] 743 744
.. ..$ : int [1:2] 1375 1376
.. ..$ : int [1:2] 335 336
.. ..$ : int [1:2] 903 904
.. ..$ : int [1:2] 651 652
.. ..$ : int [1:2] 1061 1062
.. ..$ : int [1:2] 1219 1220
.. ..$ : int [1:2] 493 494
.. ..$ : int [1:2] 745 746
.. ..$ : int [1:2] 1377 1378
.. ..$ : int [1:2] 337 338
.. ..$ : int [1:2] 905 906
.. ..$ : int [1:2] 653 654
.. ..$ : int [1:2] 1063 1064
.. ..$ : int [1:2] 1221 1222
.. ..$ : int [1:2] 495 496
.. ..$ : int [1:2] 747 748
.. ..$ : int [1:2] 1379 1380
.. ..$ : int [1:2] 339 340
.. ..$ : int [1:2] 907 908
.. ..$ : int [1:2] 655 656
.. ..$ : int [1:2] 1065 1066
.. ..$ : int [1:2] 1223 1224
.. ..$ : int [1:2] 497 498
.. ..$ : int [1:2] 749 750
.. ..$ : int [1:2] 1381 1382
.. ..$ : int [1:2] 341 342
.. ..$ : int [1:2] 909 910
.. ..$ : int [1:2] 657 658
.. .. [list output truncated]
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
# Group by Sample and Target, then summarize to get unique rows for each sample
grouped_df <- combined.fitered_df%>%
group_by(Sample, Target) %>%
summarize(Cq.Mean = mean(Cq, na.rm = TRUE)) %>%
ungroup()
str(grouped_df)
tibble [760 × 3] (S3: tbl_df/tbl/data.frame)
$ Sample : chr [1:760] "200" "200" "200" "200" ...
$ Target : chr [1:760] "ATPsynthase" "DNMT1" "GAPDH" "HSP70" ...
$ Cq.Mean: num [1:760] 25.2 33.6 25.4 31.8 26 ...
# Initialize new columns
grouped_df <- grouped_df %>%
mutate(life.stage = NA_character_,
conditioning.treatment = NA_character_,
acute.treatment = NA_character_)
# Loop through each vector
for (vec_name in names(groups_list)) {
vec <- groups_list[[vec_name]]
stage <- strsplit(vec_name, "\\.")[[1]][1]
conditioning_treatment <- strsplit(vec_name, "\\.")[[1]][2]
acute_treatment <- strsplit(vec_name, "\\.")[[1]][3]
# Loop through each row in grouped_df
for (i in 1:nrow(grouped_df)) {
sample <- grouped_df$Sample[i]
# Check if sample is in the vector
if (sample %in% vec) {
# Update life.stage and treatment columns
grouped_df$life.stage[i] <- stage
grouped_df$conditioning.treatment[i] <- conditioning_treatment
grouped_df$acute.treatment[i] <-acute_treatment
}
}
}
str(grouped_df)
tibble [760 × 6] (S3: tbl_df/tbl/data.frame)
$ Sample : chr [1:760] "200" "200" "200" "200" ...
$ Target : chr [1:760] "ATPsynthase" "DNMT1" "GAPDH" "HSP70" ...
$ Cq.Mean : num [1:760] 25.2 33.6 25.4 31.8 26 ...
$ life.stage : chr [1:760] "seed" "seed" "seed" "seed" ...
$ conditioning.treatment: chr [1:760] "treated" "treated" "treated" "treated" ...
$ acute.treatment : chr [1:760] "high" "high" "high" "high" ...
# Calculate delta Cq by subtracting GAPDH Cq.Mean from each corresponding Sample Cq.Mean
delta_Cq_df <- calculate_delta_Cq(grouped_df)
# Filters out normalizing gene, since no need to compare normalizing gene to itself.
delta_Cq_df <- delta_Cq_df %>%
filter(!is.na(life.stage), !is.na(Target), Target != "GAPDH")
str(delta_Cq_df)
tibble [665 × 7] (S3: tbl_df/tbl/data.frame)
$ Sample : chr [1:665] "200" "200" "200" "200" ...
$ Target : chr [1:665] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
$ Cq.Mean : num [1:665] 25.2 33.6 31.8 26 30.6 ...
$ life.stage : chr [1:665] "seed" "seed" "seed" "seed" ...
$ conditioning.treatment: chr [1:665] "treated" "treated" "treated" "treated" ...
$ acute.treatment : chr [1:665] "high" "high" "high" "high" ...
$ delta_Cq : num [1:665] -0.243 8.168 6.349 0.578 5.135 ...
This code does the following:
# Extract unique life.stage levels
unique_life_stages <- unique(delta_Cq_df$life.stage)
# Generate all possible pairs of life.stage levels
life_stage_pairs <- combn(unique_life_stages, 2, simplify = FALSE)
# Initialize a list to store results
life_stage_t_test_results_list <- list()
for (pair in life_stage_pairs) {
stage1 <- pair[1]
stage2 <- pair[2]
# Perform t-test for each Target comparing the two life.stage levels
t_test_results <- delta_Cq_df %>%
filter(life.stage %in% c(stage1, stage2)) %>%
group_by(Target) %>%
summarise(
t_test_result = list(t.test(delta_Cq ~ life.stage))
) %>%
ungroup() %>%
mutate(
estimate_diff = sapply(t_test_result, function(x) x$estimate[1] - x$estimate[2]),
p_value = sapply(t_test_result, function(x) x$p.value),
asterisk = ifelse(p_value <= 0.05, "*", ""), # Adds asterisk column and asterisk for p-value.
comparison = paste(stage1, "vs", stage2, sep = ".")
) %>%
select(!t_test_result)
life_stage_t_test_results_list[[paste(stage1, stage2, sep = ".")]] <- t_test_results
}
# Combine results into a single data frame
life_stage_t_test_results_df <- bind_rows(life_stage_t_test_results_list, .id = "comparison")
# View the results
print(life_stage_t_test_results_df)
# A tibble: 42 × 5
Target estimate_diff p_value asterisk comparison
<chr> <dbl> <dbl> <chr> <chr>
1 ATPsynthase 0.310 0.0183 "*" seed.adult
2 DNMT1 -0.0736 0.829 "" seed.adult
3 HSP70 0.179 0.702 "" seed.adult
4 HSP90 0.727 0.00635 "*" seed.adult
5 cGAS -0.0136 0.963 "" seed.adult
6 citrate synthase -0.215 0.391 "" seed.adult
7 viperin -1.08 0.0000737 "*" seed.adult
8 ATPsynthase 0.207 0.0819 "" seed.juvenile
9 DNMT1 -0.434 0.254 "" seed.juvenile
10 HSP70 0.657 0.179 "" seed.juvenile
# ℹ 32 more rows
This code does the following:
# Extract unique life.stage levels
unique_life_stages <- unique(delta_Cq_df$life.stage)
# Initialize a list to store results
conditioning_treatment_t_test_results_list <- list()
for (stage in unique_life_stages) {
# Extract unique conditioning.treatment levels within the current life.stage
unique_treatments <- unique(delta_Cq_df %>% filter(life.stage == stage) %>% pull(conditioning.treatment))
# Generate all possible pairs of conditioning.treatment levels
treatment_pairs <- combn(unique_treatments, 2, simplify = FALSE)
for (pair in treatment_pairs) {
treatment1 <- pair[1]
treatment2 <- pair[2]
# Perform t-test for each Target comparing the two conditioning.treatment levels within the current life.stage
t_test_results <- delta_Cq_df %>%
filter(life.stage == stage, conditioning.treatment %in% c(treatment1, treatment2)) %>%
group_by(Target) %>%
summarise(
t_test_result = list(t.test(delta_Cq ~ conditioning.treatment))
) %>%
ungroup() %>%
mutate(
estimate_diff = sapply(t_test_result, function(x) x$estimate[1] - x$estimate[2]),
p_value = sapply(t_test_result, function(x) x$p.value),
asterisk = ifelse(p_value <= 0.05, "*", ""), # Adds asterisk column and asterisk for p-value.
comparison = paste(stage, treatment1, "vs", treatment2, sep = ".")
) %>%
select(!t_test_result)
conditioning_treatment_t_test_results_list[[paste(stage, treatment1, treatment2, sep = ".")]] <- t_test_results
}
}
# Combine results into a single data frame
conditioning_treatment_t_test_results_df <- bind_rows(conditioning_treatment_t_test_results_list, .id = "comparison")
# View the results
print(conditioning_treatment_t_test_results_df)
# A tibble: 28 × 5
Target estimate_diff p_value asterisk comparison
<chr> <dbl> <dbl> <chr> <chr>
1 ATPsynthase -0.0248 0.907 "" seed.treated.control
2 DNMT1 0.206 0.746 "" seed.treated.control
3 HSP70 -0.123 0.851 "" seed.treated.control
4 HSP90 -0.299 0.439 "" seed.treated.control
5 cGAS -0.0776 0.888 "" seed.treated.control
6 citrate synthase -0.0148 0.976 "" seed.treated.control
7 viperin 0.160 0.695 "" seed.treated.control
8 ATPsynthase 0.0779 0.603 "" adult.treated.control
9 DNMT1 0.312 0.278 "" adult.treated.control
10 HSP70 -0.941 0.177 "" adult.treated.control
# ℹ 18 more rows
This code does the following:
Excludes seed and spat, as these were only held at ambient for the acute treatment.
# Extract unique life.stage levels, excluding 'seed' and 'spat'
unique_life_stages <- unique(delta_Cq_df$life.stage)
unique_life_stages <- setdiff(unique_life_stages, c("seed", "spat"))
# Initialize a list to store results
acute_treatment_t_test_results_list <- list()
for (stage in unique_life_stages) {
# Extract unique acute.treatment levels within the current life.stage
unique_treatments <- unique(delta_Cq_df %>% filter(life.stage == stage) %>% pull(acute.treatment))
# Check if there are at least 2 unique treatments
if (length(unique_treatments) >= 2) {
# Generate all possible pairs of acute.treatment levels
treatment_pairs <- combn(unique_treatments, 2, simplify = FALSE)
for (pair in treatment_pairs) {
treatment1 <- pair[1]
treatment2 <- pair[2]
# Perform t-test for each Target comparing the two acute.treatment levels within the current life.stage
t_test_results <- delta_Cq_df %>%
filter(life.stage == stage, acute.treatment %in% c(treatment1, treatment2)) %>%
group_by(Target) %>%
summarise(
t_test_result = list(t.test(delta_Cq ~ acute.treatment))
) %>%
ungroup() %>%
mutate(
estimate_diff = sapply(t_test_result, function(x) x$estimate[1] - x$estimate[2]),
p_value = sapply(t_test_result, function(x) x$p.value),
asterisk = ifelse(p_value <= 0.05, "*", ""), # Adds asterisk column and asterisk for p-value.
comparison = paste(stage, treatment1, "vs", treatment2, sep = ".")
) %>%
select(!t_test_result)
acute_treatment_t_test_results_list[[paste(stage, treatment1, treatment2, sep = ".")]] <- t_test_results
}
}
}
# Combine results into a single data frame
acute_treatment_t_test_results_df <- bind_rows(acute_treatment_t_test_results_list, .id = "comparison")
# View the results
print(acute_treatment_t_test_results_df)
# A tibble: 14 × 5
Target estimate_diff p_value asterisk comparison
<chr> <dbl> <dbl> <chr> <chr>
1 ATPsynthase 0.0605 0.687 "" adult.ambient.high
2 DNMT1 0.314 0.275 "" adult.ambient.high
3 HSP70 0.276 0.696 "" adult.ambient.high
4 HSP90 0.499 0.149 "" adult.ambient.high
5 cGAS 0.329 0.200 "" adult.ambient.high
6 citrate synthase 0.0668 0.622 "" adult.ambient.high
7 viperin 0.323 0.251 "" adult.ambient.high
8 ATPsynthase -0.0370 0.738 "" juvenile.ambient.high
9 DNMT1 -0.672 0.121 "" juvenile.ambient.high
10 HSP70 0.745 0.319 "" juvenile.ambient.high
11 HSP90 0.0450 0.859 "" juvenile.ambient.high
12 cGAS -0.203 0.474 "" juvenile.ambient.high
13 citrate synthase 0.0399 0.870 "" juvenile.ambient.high
14 viperin -0.424 0.304 "" juvenile.ambient.high
# Extract unique life.stage levels, excluding 'seed' and 'spat'
unique_life_stages <- unique(delta_Cq_df$life.stage)
#unique_life_stages <- setdiff(unique_life_stages, c("seed", "spat"))
# Extract unique conditioning.treatment levels
unique_conditioning_treatments <- unique(delta_Cq_df$conditioning.treatment)
# Initialize a list to store results
acute_treatment_within_life.stages_conditioning_t_test_results_list <- list()
for (stage in unique_life_stages) {
for (conditioning in unique_conditioning_treatments) {
# Extract unique acute.treatment levels within the current life.stage and conditioning.treatment
unique_treatments <- unique(delta_Cq_df %>% filter(life.stage == stage, conditioning.treatment == conditioning) %>% pull(acute.treatment))
# Check if there are at least 2 unique treatments
if (length(unique_treatments) >= 2) {
# Generate all possible pairs of acute.treatment levels
treatment_pairs <- combn(unique_treatments, 2, simplify = FALSE)
for (pair in treatment_pairs) {
treatment1 <- pair[1]
treatment2 <- pair[2]
# Perform t-test for each Target comparing the two acute.treatment levels within the current life.stage and conditioning.treatment
t_test_results <- delta_Cq_df %>%
filter(life.stage == stage, conditioning.treatment == conditioning, acute.treatment %in% c(treatment1, treatment2)) %>%
group_by(Target) %>%
summarise(
t_test_result = list(t.test(delta_Cq ~ acute.treatment))
) %>%
ungroup() %>%
mutate(
estimate_diff = sapply(t_test_result, function(x) x$estimate[1] - x$estimate[2]),
p_value = sapply(t_test_result, function(x) x$p.value),
asterisk = ifelse(p_value <= 0.05, "*", ""), # Adds asterisk column and asterisk for p-value.
comparison = paste(stage, conditioning, treatment1, "vs", treatment2, sep = ".")
) %>%
select(!t_test_result)
acute_treatment_within_life.stages_conditioning_t_test_results_list[[paste(stage, conditioning, treatment1, treatment2, sep = ".")]] <- t_test_results
}
}
}
}
# Combine results into a single data frame
acute_treatment_within_life.stages_conditioning_t_test_results_df <- bind_rows(acute_treatment_within_life.stages_conditioning_t_test_results_list, .id = "comparison_id")
# View the results
print(acute_treatment_within_life.stages_conditioning_t_test_results_df)
# A tibble: 56 × 6
comparison_id Target estimate_diff p_value asterisk comparison
<chr> <chr> <dbl> <dbl> <chr> <chr>
1 seed.treated.high.ambient ATPsynth… 0.658 0.0602 "" seed.trea…
2 seed.treated.high.ambient DNMT1 -1.55 0.166 "" seed.trea…
3 seed.treated.high.ambient HSP70 -0.133 0.920 "" seed.trea…
4 seed.treated.high.ambient HSP90 0.970 0.0699 "" seed.trea…
5 seed.treated.high.ambient cGAS -0.825 0.0627 "" seed.trea…
6 seed.treated.high.ambient citrate … -1.15 0.0114 "*" seed.trea…
7 seed.treated.high.ambient viperin -0.771 0.255 "" seed.trea…
8 seed.control.ambient.high ATPsynth… -0.304 0.125 "" seed.cont…
9 seed.control.ambient.high DNMT1 -1.54 0.279 "" seed.cont…
10 seed.control.ambient.high HSP70 -1.59 0.0190 "*" seed.cont…
# ℹ 46 more rows
# Create box plots for each comparison
unique_comparisons <- unique(life_stage_t_test_results_df$comparison)
for (comparison in unique_comparisons) {
create_boxplot_delta_Cq(delta_Cq_df, comparison, life_stage_t_test_results_df)
}
# Create box plots for each comparison
unique_comparisons <- unique(conditioning_treatment_t_test_results_df$comparison)
for (comparison in unique_comparisons) {
create_boxplot_conditioning(delta_Cq_df, comparison, conditioning_treatment_t_test_results_df)
}
# Create box plots for each comparison
unique_comparisons <- unique(acute_treatment_t_test_results_df$comparison)
for (comparison in unique_comparisons) {
create_boxplot_acute(delta_Cq_df, comparison, acute_treatment_t_test_results_df)
}
# Loop through each comparison in the t-test results and create box plots
for (comparison in unique(acute_treatment_within_life.stages_conditioning_t_test_results_df$comparison)) {
create_boxplot_acute_conditioning(delta_Cq_df, comparison, acute_treatment_within_life.stages_conditioning_t_test_results_df)
}
# Calculate delta_delta_Cq
delta_delta_conditioning_fold_change <- delta_Cq_df %>%
group_by(life.stage, Target) %>%
summarize(
treated_delta_Cq = mean(delta_Cq[conditioning.treatment == "treated"], na.rm = TRUE),
control_delta_Cq = mean(delta_Cq[conditioning.treatment == "control"], na.rm = TRUE)
) %>%
mutate(delta_delta_Cq = treated_delta_Cq - control_delta_Cq) %>%
select(life.stage, Target, delta_delta_Cq)
str(delta_delta_conditioning_fold_change)
gropd_df [28 × 3] (S3: grouped_df/tbl_df/tbl/data.frame)
$ life.stage : chr [1:28] "adult" "adult" "adult" "adult" ...
$ Target : chr [1:28] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
$ delta_delta_Cq: num [1:28] -0.0779 -0.3116 0.941 0.7639 0.1955 ...
- attr(*, "groups")= tibble [4 × 2] (S3: tbl_df/tbl/data.frame)
..$ life.stage: chr [1:4] "adult" "juvenile" "seed" "spat"
..$ .rows : list<int> [1:4]
.. ..$ : int [1:7] 1 2 3 4 5 6 7
.. ..$ : int [1:7] 8 9 10 11 12 13 14
.. ..$ : int [1:7] 15 16 17 18 19 20 21
.. ..$ : int [1:7] 22 23 24 25 26 27 28
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
# Calculate delta_delta_Cq for acute treatment
delta_delta_Cq_acute_df <- delta_Cq_df %>%
group_by(life.stage, Target, acute.treatment) %>%
summarize(
treated_delta_Cq = mean(delta_Cq[conditioning.treatment == "treated"], na.rm = TRUE),
control_delta_Cq = mean(delta_Cq[conditioning.treatment == "control"], na.rm = TRUE)
) %>%
mutate(delta_delta_Cq = treated_delta_Cq - control_delta_Cq) %>%
select(life.stage, Target, acute.treatment, delta_delta_Cq)
str(delta_delta_Cq_acute_df)
gropd_df [56 × 4] (S3: grouped_df/tbl_df/tbl/data.frame)
$ life.stage : chr [1:56] "adult" "adult" "adult" "adult" ...
$ Target : chr [1:56] "ATPsynthase" "ATPsynthase" "DNMT1" "DNMT1" ...
$ acute.treatment: chr [1:56] "ambient" "high" "ambient" "high" ...
$ delta_delta_Cq : num [1:56] -0.112 -0.0438 -0.2467 -0.3765 0.9455 ...
- attr(*, "groups")= tibble [28 × 3] (S3: tbl_df/tbl/data.frame)
..$ life.stage: chr [1:28] "adult" "adult" "adult" "adult" ...
..$ Target : chr [1:28] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
..$ .rows : list<int> [1:28]
.. ..$ : int [1:2] 1 2
.. ..$ : int [1:2] 3 4
.. ..$ : int [1:2] 5 6
.. ..$ : int [1:2] 7 8
.. ..$ : int [1:2] 9 10
.. ..$ : int [1:2] 11 12
.. ..$ : int [1:2] 13 14
.. ..$ : int [1:2] 15 16
.. ..$ : int [1:2] 17 18
.. ..$ : int [1:2] 19 20
.. ..$ : int [1:2] 21 22
.. ..$ : int [1:2] 23 24
.. ..$ : int [1:2] 25 26
.. ..$ : int [1:2] 27 28
.. ..$ : int [1:2] 29 30
.. ..$ : int [1:2] 31 32
.. ..$ : int [1:2] 33 34
.. ..$ : int [1:2] 35 36
.. ..$ : int [1:2] 37 38
.. ..$ : int [1:2] 39 40
.. ..$ : int [1:2] 41 42
.. ..$ : int [1:2] 43 44
.. ..$ : int [1:2] 45 46
.. ..$ : int [1:2] 47 48
.. ..$ : int [1:2] 49 50
.. ..$ : int [1:2] 51 52
.. ..$ : int [1:2] 53 54
.. ..$ : int [1:2] 55 56
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
# Calculate delta_delta_Cq for life stage comparisons
delta_delta_Cq_life_stage_df <- delta_Cq_df %>%
group_by(Target, life.stage) %>%
summarize(mean_delta_Cq = mean(delta_Cq, na.rm = TRUE)) %>%
ungroup() %>%
pivot_wider(names_from = life.stage, values_from = mean_delta_Cq) %>%
mutate(
delta_delta_Cq_adult_vs_seed = adult - seed,
delta_delta_Cq_spat_vs_seed = spat - seed,
delta_delta_Cq_adult_vs_spat = adult - spat
) %>%
pivot_longer(cols = starts_with("delta_delta_Cq_"), names_to = "comparison", values_to = "delta_delta_Cq") %>%
filter(!is.na(delta_delta_Cq))
# Display the structure of the resulting data frame
str(delta_delta_Cq_life_stage_df)
tibble [21 × 7] (S3: tbl_df/tbl/data.frame)
$ Target : chr [1:21] "ATPsynthase" "ATPsynthase" "ATPsynthase" "DNMT1" ...
$ adult : num [1:21] 0.48 0.48 0.48 6.17 6.17 ...
$ juvenile : num [1:21] 0.378 0.378 0.378 5.808 5.808 ...
$ seed : num [1:21] 0.17 0.17 0.17 6.24 6.24 ...
$ spat : num [1:21] 0.337 0.337 0.337 6.443 6.443 ...
$ comparison : chr [1:21] "delta_delta_Cq_adult_vs_seed" "delta_delta_Cq_spat_vs_seed" "delta_delta_Cq_adult_vs_spat" "delta_delta_Cq_adult_vs_seed" ...
$ delta_delta_Cq: num [1:21] 0.3098 0.1664 0.1433 -0.0736 0.2019 ...
# Calculate delta_delta_Cq for acute treatment comparisons within each life stage and conditioning treatment
delta_delta_Cq_acute_within_life_stage_conditioning_df <- delta_Cq_df %>%
group_by(life.stage, conditioning.treatment, Target, acute.treatment) %>%
summarize(mean_delta_Cq = mean(delta_Cq, na.rm = TRUE)) %>%
ungroup() %>%
pivot_wider(names_from = acute.treatment, values_from = mean_delta_Cq) %>%
mutate(delta_delta_Cq_high_vs_ambient = high - ambient) %>%
pivot_longer(cols = starts_with("delta_delta_Cq_"), names_to = "comparison", values_to = "delta_delta_Cq") %>%
filter(!is.na(delta_delta_Cq))
# Display the structure of the resulting data frame
str(delta_delta_Cq_acute_within_life_stage_conditioning_df)
tibble [56 × 7] (S3: tbl_df/tbl/data.frame)
$ life.stage : chr [1:56] "adult" "adult" "adult" "adult" ...
$ conditioning.treatment: chr [1:56] "control" "control" "control" "control" ...
$ Target : chr [1:56] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
$ ambient : num [1:56] 0.566 6.448 3.944 1.259 5.207 ...
$ high : num [1:56] 0.472 6.199 3.673 0.29 4.609 ...
$ comparison : chr [1:56] "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" ...
$ delta_delta_Cq : num [1:56] -0.0946 -0.2492 -0.2715 -0.969 -0.5983 ...
# Calculate fold change and output to a new data frame
fold_change_life_stage_df <- delta_delta_Cq_life_stage_df %>%
mutate(fold_change = 2^(-delta_delta_Cq))
# Display the structure of the resulting data frame
str(fold_change_life_stage_df)
tibble [21 × 8] (S3: tbl_df/tbl/data.frame)
$ Target : chr [1:21] "ATPsynthase" "ATPsynthase" "ATPsynthase" "DNMT1" ...
$ adult : num [1:21] 0.48 0.48 0.48 6.17 6.17 ...
$ juvenile : num [1:21] 0.378 0.378 0.378 5.808 5.808 ...
$ seed : num [1:21] 0.17 0.17 0.17 6.24 6.24 ...
$ spat : num [1:21] 0.337 0.337 0.337 6.443 6.443 ...
$ comparison : chr [1:21] "delta_delta_Cq_adult_vs_seed" "delta_delta_Cq_spat_vs_seed" "delta_delta_Cq_adult_vs_spat" "delta_delta_Cq_adult_vs_seed" ...
$ delta_delta_Cq: num [1:21] 0.3098 0.1664 0.1433 -0.0736 0.2019 ...
$ fold_change : num [1:21] 0.807 0.891 0.905 1.052 0.869 ...
delta_delta_conditioning_fold_change <- delta_delta_conditioning_fold_change %>%
mutate(fold_change = 2^(-delta_delta_Cq)) %>%
distinct(Target, fold_change)
str(delta_delta_conditioning_fold_change)
gropd_df [28 × 3] (S3: grouped_df/tbl_df/tbl/data.frame)
$ life.stage : chr [1:28] "adult" "adult" "adult" "adult" ...
$ Target : chr [1:28] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
$ fold_change: num [1:28] 1.055 1.241 0.521 0.589 0.873 ...
- attr(*, "groups")= tibble [4 × 2] (S3: tbl_df/tbl/data.frame)
..$ life.stage: chr [1:4] "adult" "juvenile" "seed" "spat"
..$ .rows : list<int> [1:4]
.. ..$ : int [1:7] 1 2 3 4 5 6 7
.. ..$ : int [1:7] 8 9 10 11 12 13 14
.. ..$ : int [1:7] 15 16 17 18 19 20 21
.. ..$ : int [1:7] 22 23 24 25 26 27 28
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
# Calculate fold change for acute treatment
delta_delta_acute_fold_change <- delta_delta_Cq_acute_df %>%
mutate(fold_change = 2^(-delta_delta_Cq)) %>%
distinct(life.stage, Target, acute.treatment, fold_change)
# Display the structure of the resulting data frame
str(delta_delta_acute_fold_change)
gropd_df [56 × 4] (S3: grouped_df/tbl_df/tbl/data.frame)
$ life.stage : chr [1:56] "adult" "adult" "adult" "adult" ...
$ Target : chr [1:56] "ATPsynthase" "ATPsynthase" "DNMT1" "DNMT1" ...
$ acute.treatment: chr [1:56] "ambient" "high" "ambient" "high" ...
$ fold_change : num [1:56] 1.081 1.031 1.186 1.298 0.519 ...
- attr(*, "groups")= tibble [28 × 3] (S3: tbl_df/tbl/data.frame)
..$ life.stage: chr [1:28] "adult" "adult" "adult" "adult" ...
..$ Target : chr [1:28] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
..$ .rows : list<int> [1:28]
.. ..$ : int [1:2] 1 2
.. ..$ : int [1:2] 3 4
.. ..$ : int [1:2] 5 6
.. ..$ : int [1:2] 7 8
.. ..$ : int [1:2] 9 10
.. ..$ : int [1:2] 11 12
.. ..$ : int [1:2] 13 14
.. ..$ : int [1:2] 15 16
.. ..$ : int [1:2] 17 18
.. ..$ : int [1:2] 19 20
.. ..$ : int [1:2] 21 22
.. ..$ : int [1:2] 23 24
.. ..$ : int [1:2] 25 26
.. ..$ : int [1:2] 27 28
.. ..$ : int [1:2] 29 30
.. ..$ : int [1:2] 31 32
.. ..$ : int [1:2] 33 34
.. ..$ : int [1:2] 35 36
.. ..$ : int [1:2] 37 38
.. ..$ : int [1:2] 39 40
.. ..$ : int [1:2] 41 42
.. ..$ : int [1:2] 43 44
.. ..$ : int [1:2] 45 46
.. ..$ : int [1:2] 47 48
.. ..$ : int [1:2] 49 50
.. ..$ : int [1:2] 51 52
.. ..$ : int [1:2] 53 54
.. ..$ : int [1:2] 55 56
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
# Calculate fold change for acute treatment comparisons within each life stage and conditioning treatment
fold_change_acute_within_life_stage_conditioning_df <- delta_delta_Cq_acute_within_life_stage_conditioning_df %>%
mutate(fold_change = 2^(-delta_delta_Cq))
# Display the structure of the resulting data frame
str(fold_change_acute_within_life_stage_conditioning_df)
tibble [56 × 8] (S3: tbl_df/tbl/data.frame)
$ life.stage : chr [1:56] "adult" "adult" "adult" "adult" ...
$ conditioning.treatment: chr [1:56] "control" "control" "control" "control" ...
$ Target : chr [1:56] "ATPsynthase" "DNMT1" "HSP70" "HSP90" ...
$ ambient : num [1:56] 0.566 6.448 3.944 1.259 5.207 ...
$ high : num [1:56] 0.472 6.199 3.673 0.29 4.609 ...
$ comparison : chr [1:56] "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" "delta_delta_Cq_high_vs_ambient" ...
$ delta_delta_Cq : num [1:56] -0.0946 -0.2492 -0.2715 -0.969 -0.5983 ...
$ fold_change : num [1:56] 1.07 1.19 1.21 1.96 1.51 ...
library(ggplot2)
# Generate bar plots for each group of comparison within each life stage and conditioning treatment
plot_list <- fold_change_acute_within_life_stage_conditioning_df %>%
split(list(.$life.stage, .$conditioning.treatment, .$comparison)) %>%
lapply(function(df) {
life_stage <- unique(df$life.stage)
conditioning_treatment <- unique(df$conditioning.treatment)
comparison_title <- gsub("delta_delta_Cq_", "", unique(df$comparison))
comparison_title <- gsub("_vs_", " vs. ", comparison_title)
ggplot(df, aes(x = Target, y = fold_change)) +
geom_bar(stat = "identity") +
labs(title = paste("Gene Expression -", life_stage, "-", conditioning_treatment, "-", comparison_title),
x = "Target", y = "Fold Change") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
})
# Display the plots
for (plot in plot_list) {
print(plot)
}
# Generate bar plots for each group of comparison
plot_list <- fold_change_life_stage_df %>%
split(.$comparison) %>%
lapply(function(df) {
comparison_title <- gsub("delta_delta_Cq_", "", unique(df$comparison))
comparison_title <- gsub("_vs_", " vs. ", comparison_title)
ggplot(df, aes(x = Target, y = fold_change)) +
geom_bar(stat = "identity") +
labs(title = paste("Gene Expression -", comparison_title), x = "Target", y = "Fold Change") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
})
# Display the plots
for (plot in plot_list) {
print(plot)
}
hist(delta_Cq_df$delta_Cq)
Run an anova model to test for effects of lifestage, conditioning, and acute treatment on delta Cq values for each target.
ATP synthase is an enzyme complex that functions to synthesize adenosine triphosphate (ATP) from adenosine diphosphate (ADP) and inorganic phosphate (Pi), essentially generating the cell’s primary energy currency by harnessing the energy from a proton gradient across a membrane.
library(car)
library(emmeans)
model<-delta_Cq_df%>%
filter(Target=="ATPsynthase")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
Df Sum Sq Mean Sq F value
life.stage 3 1.127 0.3755 2.589
conditioning.treatment 1 0.118 0.1180 0.814
acute.treatment 1 0.367 0.3665 2.527
life.stage:conditioning.treatment 3 0.116 0.0387 0.267
life.stage:acute.treatment 3 0.487 0.1624 1.119
conditioning.treatment:acute.treatment 1 0.264 0.2640 1.820
life.stage:conditioning.treatment:acute.treatment 3 0.745 0.2482 1.711
Residuals 79 11.460 0.1451
Pr(>F)
life.stage 0.0588 .
conditioning.treatment 0.3698
acute.treatment 0.1159
life.stage:conditioning.treatment 0.8492
life.stage:acute.treatment 0.3464
conditioning.treatment:acute.treatment 0.1812
life.stage:conditioning.treatment:acute.treatment 0.1715
Residuals
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(model$residuals)
[1] 31 77
leveneTest(model$residuals ~ life.stage*conditioning.treatment*acute.treatment, data=delta_Cq_df%>%filter(Target=="ATPsynthase"))
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 15 0.7829 0.6923
79
No significant effects.
The DNMT1 gene provides instructions for making an enzyme called DNA methyltransferase 1. This enzyme is involved in DNA methylation, which is the addition of methyl groups, consisting of one carbon atom and three hydrogen atoms, to DNA molecules.
model<-delta_Cq_df%>%
filter(Target=="DNMT1")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
Df Sum Sq Mean Sq F value
life.stage 3 4.78 1.592 1.296
conditioning.treatment 1 0.20 0.201 0.164
acute.treatment 1 3.87 3.869 3.149
life.stage:conditioning.treatment 3 6.67 2.222 1.809
life.stage:acute.treatment 3 10.00 3.334 2.714
conditioning.treatment:acute.treatment 1 0.96 0.959 0.781
life.stage:conditioning.treatment:acute.treatment 3 6.80 2.266 1.845
Residuals 79 97.06 1.229
Pr(>F)
life.stage 0.2817
conditioning.treatment 0.6868
acute.treatment 0.0798 .
life.stage:conditioning.treatment 0.1524
life.stage:acute.treatment 0.0504 .
conditioning.treatment:acute.treatment 0.3796
life.stage:conditioning.treatment:acute.treatment 0.1459
Residuals
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(model$residuals)
[1] 7 21
leveneTest(model$residuals ~ life.stage*conditioning.treatment*acute.treatment, data=delta_Cq_df%>%filter(Target=="DNMT1"))
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 15 0.8095 0.664
79
No significant effects.
Heat Shock Protein 70 (Hsp70) is a molecular chaperone that plays crucial roles in maintaining cellular protein homeostasis and protecting cells from stress.
model<-delta_Cq_df%>%
filter(Target=="HSP70")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
Df Sum Sq Mean Sq F value
life.stage 3 13.48 4.494 1.554
conditioning.treatment 1 0.45 0.445 0.154
acute.treatment 1 1.71 1.710 0.591
life.stage:conditioning.treatment 3 10.31 3.438 1.189
life.stage:acute.treatment 3 6.61 2.203 0.762
conditioning.treatment:acute.treatment 1 21.29 21.286 7.361
life.stage:conditioning.treatment:acute.treatment 3 16.05 5.349 1.850
Residuals 79 228.45 2.892
Pr(>F)
life.stage 0.20717
conditioning.treatment 0.69582
acute.treatment 0.44418
life.stage:conditioning.treatment 0.31938
life.stage:acute.treatment 0.51889
conditioning.treatment:acute.treatment 0.00818 **
life.stage:conditioning.treatment:acute.treatment 0.14496
Residuals
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(model$residuals)
[1] 94 9
leveneTest(model$residuals ~ life.stage*conditioning.treatment*acute.treatment, data=delta_Cq_df%>%filter(Target=="HSP70"))
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 15 0.61 0.8584
79
emm<-emmeans(model, ~ conditioning.treatment:acute.treatment | life.stage)
pairs(emm)
life.stage = adult:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -0.946 0.850 79 -1.112 0.6832
control ambient - control high 0.271 0.850 79 0.319 0.9887
control ambient - treated high -0.665 0.850 79 -0.782 0.8624
treated ambient - control high 1.217 0.850 79 1.431 0.4838
treated ambient - treated high 0.281 0.850 79 0.330 0.9875
control high - treated high -0.936 0.850 79 -1.101 0.6898
life.stage = juvenile:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -1.304 0.982 79 -1.329 0.5477
control ambient - control high -0.603 0.982 79 -0.614 0.9274
control ambient - treated high 0.788 0.982 79 0.803 0.8529
treated ambient - control high 0.702 0.982 79 0.715 0.8910
treated ambient - treated high 2.093 0.982 79 2.131 0.1522
control high - treated high 1.391 0.982 79 1.417 0.4928
life.stage = seed:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -0.608 0.982 79 -0.619 0.9257
control ambient - control high -1.589 1.202 79 -1.321 0.5523
control ambient - treated high -0.741 1.202 79 -0.616 0.9267
treated ambient - control high -0.981 1.202 79 -0.816 0.8469
treated ambient - treated high -0.133 1.202 79 -0.110 0.9995
control high - treated high 0.848 1.388 79 0.611 0.9284
life.stage = spat:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -1.236 1.030 79 -1.200 0.6286
control ambient - control high -1.873 1.098 79 -1.707 0.3271
control ambient - treated high 1.255 0.982 79 1.278 0.5794
treated ambient - control high -0.638 1.141 79 -0.559 0.9438
treated ambient - treated high 2.491 1.030 79 2.419 0.0817
control high - treated high 3.128 1.098 79 2.850 0.0280
P value adjustment: tukey method for comparing a family of 4 estimates
Significant effect of conditioning x acute treatment.
Heat shock protein 90 (Hsp90) is a molecular chaperone that helps proteins fold, mature, and remain active. Hsp90 also helps regulate signaling networks and is involved in many cellular processes.
model<-delta_Cq_df%>%
filter(Target=="HSP90")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
Df Sum Sq Mean Sq F value
life.stage 3 22.10 7.368 14.629
conditioning.treatment 1 0.98 0.983 1.952
acute.treatment 1 12.31 12.308 24.438
life.stage:conditioning.treatment 3 4.77 1.591 3.159
life.stage:acute.treatment 3 6.58 2.192 4.353
conditioning.treatment:acute.treatment 1 0.10 0.103 0.204
life.stage:conditioning.treatment:acute.treatment 3 3.20 1.067 2.118
Residuals 79 39.79 0.504
Pr(>F)
life.stage 1.15e-07 ***
conditioning.treatment 0.16627
acute.treatment 4.22e-06 ***
life.stage:conditioning.treatment 0.02921 *
life.stage:acute.treatment 0.00686 **
conditioning.treatment:acute.treatment 0.65283
life.stage:conditioning.treatment:acute.treatment 0.10451
Residuals
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(model$residuals)
[1] 94 80
leveneTest(model$residuals ~ life.stage*conditioning.treatment*acute.treatment, data=delta_Cq_df%>%filter(Target=="HSP90"))
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 15 1.5079 0.1226
79
emm<-emmeans(model, ~ conditioning.treatment:acute.treatment | life.stage)
pairs(emm)
life.stage = adult:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -0.2939 0.355 79 -0.828 0.8409
control ambient - control high 0.9690 0.355 79 2.731 0.0382
control ambient - treated high -0.2648 0.355 79 -0.746 0.8780
treated ambient - control high 1.2630 0.355 79 3.559 0.0035
treated ambient - treated high 0.0292 0.355 79 0.082 0.9998
control high - treated high -1.2338 0.355 79 -3.477 0.0045
life.stage = juvenile:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -0.1098 0.410 79 -0.268 0.9932
control ambient - control high -0.0451 0.410 79 -0.110 0.9995
control ambient - treated high 0.0254 0.410 79 0.062 0.9999
treated ambient - control high 0.0647 0.410 79 0.158 0.9986
treated ambient - treated high 0.1352 0.410 79 0.330 0.9875
control high - treated high 0.0704 0.410 79 0.172 0.9982
life.stage = seed:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -0.1505 0.410 79 -0.367 0.9829
control ambient - control high 1.4140 0.502 79 2.818 0.0305
control ambient - treated high 0.8191 0.502 79 1.632 0.3666
treated ambient - control high 1.5645 0.502 79 3.118 0.0133
treated ambient - treated high 0.9696 0.502 79 1.932 0.2231
control high - treated high -0.5948 0.579 79 -1.027 0.7344
life.stage = spat:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -0.1549 0.430 79 -0.360 0.9839
control ambient - control high 0.8926 0.458 79 1.949 0.2165
control ambient - treated high 1.7432 0.410 79 4.255 0.0003
treated ambient - control high 1.0475 0.476 79 2.200 0.1321
treated ambient - treated high 1.8981 0.430 79 4.417 0.0002
control high - treated high 0.8506 0.458 79 1.857 0.2551
P value adjustment: tukey method for comparing a family of 4 estimates
Significant effect of lifestage x acute treatment, lifestage x conditioning treatment, acute treatment, and lifestage.
The cGAS gene is involved in several processes, including cellular response to exogenous dsRNA, positive regulation of intracellular signal transduction, and regulation of defense response.
model<-delta_Cq_df%>%
filter(Target=="cGAS")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
Df Sum Sq Mean Sq F value
life.stage 3 5.18 1.7253 2.721
conditioning.treatment 1 0.01 0.0147 0.023
acute.treatment 1 0.37 0.3707 0.585
life.stage:conditioning.treatment 3 0.57 0.1901 0.300
life.stage:acute.treatment 3 6.68 2.2251 3.510
conditioning.treatment:acute.treatment 1 0.84 0.8435 1.330
life.stage:conditioning.treatment:acute.treatment 3 4.27 1.4238 2.246
Residuals 79 50.09 0.6340
Pr(>F)
life.stage 0.0499 *
conditioning.treatment 0.8794
acute.treatment 0.4468
life.stage:conditioning.treatment 0.8254
life.stage:acute.treatment 0.0190 *
conditioning.treatment:acute.treatment 0.2522
life.stage:conditioning.treatment:acute.treatment 0.0894 .
Residuals
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(model$residuals)
[1] 70 10
leveneTest(model$residuals ~ life.stage*conditioning.treatment*acute.treatment, data=delta_Cq_df%>%filter(Target=="cGAS"))
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 15 0.817 0.656
79
emm<-emmeans(model, ~ conditioning.treatment:acute.treatment | life.stage)
pairs(emm)
life.stage = adult:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient 0.0739 0.398 79 0.186 0.9977
control ambient - control high 0.5983 0.398 79 1.503 0.4407
control ambient - treated high 0.1334 0.398 79 0.335 0.9869
treated ambient - control high 0.5244 0.398 79 1.317 0.5549
treated ambient - treated high 0.0595 0.398 79 0.150 0.9988
control high - treated high -0.4649 0.398 79 -1.168 0.6489
life.stage = juvenile:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -0.0378 0.460 79 -0.082 0.9998
control ambient - control high -0.2947 0.460 79 -0.641 0.9183
control ambient - treated high -0.1488 0.460 79 -0.324 0.9882
treated ambient - control high -0.2569 0.460 79 -0.559 0.9439
treated ambient - treated high -0.1110 0.460 79 -0.241 0.9950
control high - treated high 0.1459 0.460 79 0.317 0.9888
life.stage = seed:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -0.3407 0.460 79 -0.741 0.8802
control ambient - control high -1.6144 0.563 79 -2.867 0.0267
control ambient - treated high -1.1659 0.563 79 -2.071 0.1717
treated ambient - control high -1.2737 0.563 79 -2.262 0.1158
treated ambient - treated high -0.8252 0.563 79 -1.466 0.4629
control high - treated high 0.4485 0.650 79 0.690 0.9007
life.stage = spat:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -0.6338 0.482 79 -1.315 0.5565
control ambient - control high -0.8533 0.514 79 -1.660 0.3515
control ambient - treated high 0.2519 0.460 79 0.548 0.9468
treated ambient - control high -0.2195 0.534 79 -0.411 0.9764
treated ambient - treated high 0.8857 0.482 79 1.837 0.2639
control high - treated high 1.1053 0.514 79 2.150 0.1464
P value adjustment: tukey method for comparing a family of 4 estimates
Significant effect of lifestage x acute treatment and lifestage.
Citrate synthase is important for energy production in the TCA cycle and is linked to the electron transport chain. It is also used as an enzyme marker for intact mitochondria.
model<-delta_Cq_df%>%
filter(Target=="citrate synthase")%>%
aov(delta_Cq ~ life.stage * conditioning.treatment * acute.treatment, data=.)
summary(model)
Df Sum Sq Mean Sq F value
life.stage 3 1.383 0.461 1.985
conditioning.treatment 1 0.210 0.210 0.904
acute.treatment 1 5.025 5.025 21.641
life.stage:conditioning.treatment 3 0.146 0.049 0.210
life.stage:acute.treatment 3 10.632 3.544 15.262
conditioning.treatment:acute.treatment 1 0.012 0.012 0.052
life.stage:conditioning.treatment:acute.treatment 3 1.430 0.477 2.053
Residuals 79 18.345 0.232
Pr(>F)
life.stage 0.123
conditioning.treatment 0.345
acute.treatment 1.30e-05 ***
life.stage:conditioning.treatment 0.889
life.stage:acute.treatment 6.36e-08 ***
conditioning.treatment:acute.treatment 0.820
life.stage:conditioning.treatment:acute.treatment 0.113
Residuals
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(model$residuals)
[1] 10 93
leveneTest(model$residuals ~ life.stage*conditioning.treatment*acute.treatment, data=delta_Cq_df%>%filter(Target=="citrate synthase"))
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 15 0.7099 0.7674
79
emm<-emmeans(model, ~ conditioning.treatment:acute.treatment | life.stage)
pairs(emm)
life.stage = adult:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient 0.31754 0.241 79 1.318 0.5544
control ambient - control high 0.30936 0.241 79 1.284 0.5758
control ambient - treated high 0.14178 0.241 79 0.588 0.9353
treated ambient - control high -0.00818 0.241 79 -0.034 1.0000
treated ambient - treated high -0.17576 0.241 79 -0.729 0.8850
control high - treated high -0.16758 0.241 79 -0.696 0.8985
life.stage = juvenile:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient 0.16164 0.278 79 0.581 0.9375
control ambient - control high 0.04161 0.278 79 0.150 0.9988
control ambient - treated high 0.19990 0.278 79 0.719 0.8894
treated ambient - control high -0.12003 0.278 79 -0.431 0.9729
treated ambient - treated high 0.03826 0.278 79 0.138 0.9991
control high - treated high 0.15830 0.278 79 0.569 0.9410
life.stage = seed:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient -0.33494 0.278 79 -1.204 0.6262
control ambient - control high -2.10686 0.341 79 -6.183 <.0001
control ambient - treated high -1.48124 0.341 79 -4.347 0.0002
treated ambient - control high -1.77193 0.341 79 -5.200 <.0001
treated ambient - treated high -1.14630 0.341 79 -3.364 0.0064
control high - treated high 0.62562 0.393 79 1.590 0.3901
life.stage = spat:
contrast estimate SE df t.ratio p.value
control ambient - treated ambient 0.37761 0.292 79 1.294 0.5694
control ambient - control high -0.89118 0.311 79 -2.865 0.0269
control ambient - treated high -0.70973 0.278 79 -2.551 0.0599
treated ambient - control high -1.26879 0.323 79 -3.925 0.0010
treated ambient - treated high -1.08734 0.292 79 -3.726 0.0020
control high - treated high 0.18145 0.311 79 0.583 0.9368
P value adjustment: tukey method for comparing a family of 4 estimates
Significant effect of acute treatment and lifestage x acute treatment.
All tests pass normality and homogeneity of variance. ANOVA tests are appropriate.
Display plot of acute x conditioning treatment faceted for each lifestage for each target.
Show sample size for each group.
delta_Cq_df%>%
group_by(life.stage, acute.treatment, conditioning.treatment)%>%
summarise(n=length(unique(Sample)))
# A tibble: 16 × 4
# Groups: life.stage, acute.treatment [8]
life.stage acute.treatment conditioning.treatment n
<chr> <chr> <chr> <int>
1 adult ambient control 8
2 adult ambient treated 8
3 adult high control 8
4 adult high treated 8
5 juvenile ambient control 6
6 juvenile ambient treated 6
7 juvenile high control 6
8 juvenile high treated 6
9 seed ambient control 6
10 seed ambient treated 6
11 seed high control 3
12 seed high treated 3
13 spat ambient control 6
14 spat ambient treated 5
15 spat high control 4
16 spat high treated 6
Note that sample sizes for spat, seed, or juveniles that are <6 are due to lack of RNA isolated from the samples.
Adults were all sampled from one family (Pink). Previous flow cytometry was unclear as to diploid or triploid classification. We have extracted DNA that we can use to verify ploidy.
plot1<-delta_Cq_df%>%
filter(Target=="ATPsynthase")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_point()+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("ATP synthase (no sign. effects)")+
ylab("Delta Cq")+
ylim(-0.5,1.5)+
geom_hline(yintercept=0, linetype="dashed")+
theme_classic();plot1
plot2<-delta_Cq_df%>%
filter(Target=="DNMT1")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
geom_point()+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("DMNT1 (no sign. effects)")+
ylab("Delta Cq")+
geom_hline(yintercept=0, linetype="dashed")+
ylim(0,9)+
theme_classic();plot2
plot3<-delta_Cq_df%>%
filter(Target=="HSP70")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_point()+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("HSP70 (sign. conditioning x acute treatment)")+
ylab("Delta Cq")+
ylim(0,8)+
geom_hline(yintercept=0, linetype="dashed")+
theme_classic();plot3
plot3a<-delta_Cq_df%>%
filter(Target=="HSP70")%>%
group_by(Target, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
geom_point()+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("HSP70 (sign. conditioning x acute treatment)")+
ylab("Delta Cq")+
geom_hline(yintercept=0, linetype="dashed")+
ylim(0,8)+
theme_classic();plot3a
plot4<-delta_Cq_df%>%
filter(Target=="HSP90")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_point()+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("HSP90 (sign. lifestage x acute & lifestage x conditioning treatment)")+
ylab("Delta Cq")+
ylim(-2,2.5)+
geom_hline(yintercept=0, linetype="dashed")+
theme_classic();plot4
plot5<-delta_Cq_df%>%
filter(Target=="cGAS")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_point()+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("cGAS (sign. lifestage x acute treatment)")+
ylab("Delta Cq")+
ylim(2,8)+
geom_hline(yintercept=0, linetype="dashed")+
theme_classic();plot5
plot6<-delta_Cq_df%>%
filter(Target=="citrate synthase")%>%
group_by(Target, life.stage, acute.treatment, conditioning.treatment)%>%
summarise(mean=mean(delta_Cq, na.rm=TRUE), se=sd(delta_Cq, na.rm=TRUE)/sqrt(length(delta_Cq)))%>%
ggplot(aes(x=acute.treatment, y=mean, colour=conditioning.treatment))+
facet_grid(~life.stage)+
scale_colour_manual(values=c("darkgray", "orange"))+
geom_point()+
geom_line(aes(group=conditioning.treatment))+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
ggtitle("cGAS (sign. lifestage x acute treatment)")+
ylab("Delta Cq")+
ylim(-2,3)+
geom_hline(yintercept=0, linetype="dashed")+
theme_classic();plot6